#Alexander F. Gazmararian
#afg2@princeton.edu

library(tidyverse)
library(modelsummary)
library(lmtest)
library(sandwich)
library(kableExtra)
library(margins)

source("code/fun/book_theme.r")
source("code/fun/savefig.r")
source("code/fun/fix_txt.r")
source("code/fun/coefnames4tables.r")

# Load CivicPulse data
cp <- readRDS("data/CivicPulse_Public.rds")

# Load national sample
caps <- readRDS("data/NatCAPS_20220816.rds")
caps <- rename(caps, income = Income_long) # Rename income column

# Merge and prepare data ----
caps_biz <- subset(caps, select = c(BizPred:BizPay, sample, nweightNL33D_NEW2))
cp_biz <- subset(cp, select = c(BizPred:BizPay, Sample, Weight_1))
caps_biz <- caps_biz %>% rename(weights = nweightNL33D_NEW2)
cp_biz <- cp_biz %>% rename(weights = Weight_1, sample = Sample)
biz <- rbind(caps_biz, cp_biz)

# Prepare data for analysis
biz <-
  biz %>%
  pivot_longer(cols = c(BizPred:BizPay))
biz <-
  biz %>%
  mutate(
    sampleSlim = dplyr::case_when(
      sample %in% c("Standard", "Targeted") ~ "Policymakers",
      sample == "CAPS" ~ "National"
    ),
    sample = dplyr::case_when(
      sample == "CAPS" ~ "National",
      sample == "Standard" ~ "Policymakers (National)",
      sample == "Targeted" ~ "Policymakers (Fossil Fuel)"
    )
  )

# Figure 4.5 ----

localeff.plot <-
  biz %>%
  dplyr::filter(name != "BizPred" & !is.na(value)) %>%
  mutate(
    name = dplyr::case_when(
      name == "BizTemp" ~ "How long would the jobs last after the investment ends?",
      name == "BizLocal" ~ "How many local jobs would be created?",
      name == "BizPay" ~ "How well would the jobs pay?"
    )
  ) %>%
  mutate(value = ifelse(value %in% c("Very sure", "Somewhat sure"), 1, 0)) %>%
  group_by(name, value, sample) %>%
  summarise(n = sum(weights)) %>%
  group_by(sample, name) %>%
  mutate(pct = prop.table(n)) %>%
  group_by(sample,name) %>%
  mutate(grand_n = sum(n)) %>%
  ungroup() %>%
  mutate(
    lb = pct + qnorm(0.025) * sd(pct) / sqrt(grand_n),
    ub = pct + qnorm(0.975) * sd(pct) / sqrt(grand_n),
    lb90 = pct + qnorm(0.05) * sd(pct) / sqrt(grand_n),
    ub90 = pct + qnorm(0.95) * sd(pct) / sqrt(grand_n)
  ) %>%
  filter(value==1) %>%
  ggplot(aes(x = name, y = pct, color = sample, shape = sample, label = scales::percent(pct, accuracy = 1))) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(1)) +
  geom_errorbar(aes( ymin = lb90, ymax = ub90), width = 0, size = 1.75, position = position_dodge(1)) +
  geom_point(size = 4, position = position_dodge(1)) +
  geom_text(vjust = -1.3, size = 3, position = position_dodge(1), color = "black", aes(group = sample)) +
  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_x_discrete(labels = scales::label_wrap(20),expand = c(0,.75)) +
  labs(y = "Percent very or somewhat sure", x = "", caption = "Notes: Thin and thick bars denote 95 and 90 percent confidence intervals",
       color = "", shape = "") +
  scale_color_grey(start = 0, end = .7) +
  book_theme +
  theme(
    legend.position = "bottom",
    panel.grid.major.y  = element_line(color = "lightgrey", size = .25)
  )
localeff.plot
savefig(localeff.plot, "4.5_figure_businessbeliefs", height = 2.7, filepath = "figures/")  

# Estimate linear models for online appendix ----
## National public ----
# Local jobs
capsbizlocal <- list()
capsbizlocal[[1]] <- lm(bizlocal_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary, caps, weights = nweightNL33D_NEW2)
capsbizlocal[[2]] <- lm(bizlocal_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income, caps, weights = nweightNL33D_NEW2)
capsbizlocal[[3]] <- lm(bizlocal_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income + envt_import, caps, weights = nweightNL33D_NEW2)
capsbizlocal[[4]] <- lm(bizlocal_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income + biz_reg, caps, weights = nweightNL33D_NEW2)
# Save table
file <- "tables/ch4/ols_bizlocal_localjobs.txt"
modelsummary(capsbizlocal, vcov = "HC2", stars = c("*"=.1,"**"=.05,"***"=.01),
             gof_map = c("nobs", "adj.r.squared"),
             coef_map = coefnames,
             escape=FALSE,
             output = "latex"
) %>%
  cat(., file = file)
fix_txt(file)

# Temporary jobs 
capstemp <- list()
capstemp[[1]] <- lm(tempjobs_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary, caps, weights = nweightNL33D_NEW2)
capstemp[[2]] <- lm(tempjobs_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income, caps, weights = nweightNL33D_NEW2)
capstemp[[3]] <- lm(tempjobs_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income + envt_import, caps, weights = nweightNL33D_NEW2)
capstemp[[4]] <- lm(tempjobs_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income + biz_reg, caps, weights = nweightNL33D_NEW2)
file <- "tables/ch4/ols_bizlocal_temp.txt"
modelsummary(capstemp, vcov = "HC2", stars = c("*"=.1,"**"=.05,"***"=.01),
             gof_map = c("nobs", "adj.r.squared"),
             coef_map = coefnames,
             escape=FALSE,
             output = "latex"
) %>%
  cat(., file = file)
fix_txt(file)

# Job pay
capspay <- list()
capspay[[1]] <- lm(bizpay_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary, caps, weights = nweightNL33D_NEW2)
capspay[[2]] <- lm(bizpay_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income, caps, weights = nweightNL33D_NEW2)
capspay[[3]] <- lm(bizpay_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income + envt_import, caps, weights = nweightNL33D_NEW2)
capspay[[4]] <- lm(bizpay_scale ~ age_cat + Female + Hispanic + Black + CollegeDegree + PartySummary + income + biz_reg, caps, weights = nweightNL33D_NEW2)
file <- "tables/ch4/ols_bizlocal_pay.txt"
modelsummary(capspay, vcov = "HC2", stars = c("*"=.1,"**"=.05,"***"=.01),
             gof_map = c("nobs", "adj.r.squared"),
             coef_map = coefnames,
             escape=FALSE,
             output = "latex"
) %>%
  cat(., file = file)
fix_txt(file)

## Local policymakers ----

## NOTE: the following variables are only in the restricted access data (see README.txt): biden2020 + Census_area_college + log(Census_area_population) + Census_area_urban

f <- y ~ Age_cat + Woman + College + Minority_bin + Gov_type + biden2020 + Sample + Census_area_college + log(Census_area_population) + Census_area_urban + IdeoSimple + PartySummary
f.public <- y ~ Age_cat + Woman + College + Minority_bin + Gov_type + Sample + IdeoSimple + PartySummary

cpbiz <- list()
cpbiz[["Local"]] <- lm(update(f.public, BizLocal_scale ~ .), cp, weights = Weight_1)
cpbiz[["Temporary"]] <- lm(update(f.public, BizTemp_scale ~ .), cp, weights = Weight_1)
cpbiz[["Pay"]] <- lm(update(f.public, BizPay_scale ~ .), cp, weights = Weight_1)

# Save table
file <- "tables/ch4/ols_bizlocal_civicpulse.txt"
modelsummary(cpbiz,
             vcov = "HC2",
             stars = c("*"=.1,"**"=.05,"***"=.01),
             gof_map = c("nobs", "adj.r.squared"),
             coef_map = coefnames,
             escape=FALSE,
             output = "latex"
             ) %>%
  cat(., file = file)
fix_txt(file)

# Relative perceptions of green industries ----
relative.ols <- list()
relative.ols[[1]] <- lm(update(f.public, GreenRelative_bin ~ .), cp, weights = Weight_1)
#relative.ols[[2]] <- lm(update(f, GreenRelative_bin ~ .), cp, weights = Weight_1)
file <- "tables/ch4/ols_civicpulse_greenrelative.txt"
modelsummary(relative.ols,
             vcov = "HC2",
             stars = c("*"=.1,"**"=.05,"***"=.01),
             gof_map = c("nobs", "adj.r.squared"),
             coef_map = coefnames,
             escape=FALSE,
             output = "latex") %>%
  cat(., file = file)
fix_txt(file)
